Phase Diagram of Commensurate Two-Dimensional Disordered Bose Hubbard Model 
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We establish the full groundstate phase diagram of disordered Bose-Hubbard model in two- 
dimensions at unity filling factor via quantum Monte Carlo simulations. Similarly to the three- 
dimensional case we observe extended superfluid regions persisting up to extremely large values 
of disorder and interaction strength which, however, have small superfluid fractions and thus low 
transition temperatures. In the vicinity of the superfluid-insulator transition of the pure system, 
we observe an unexpectedly weak — almost not resolvable — sensitivity of the critical interaction to 
the strength of (weak) disorder. 
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Disordered systems keep attracting a lot of attention 
as they reveal rich and nontrivial physics of interplay be- 
tween interaction effects and localization [TJ [2] . Bosons 
are particularly interesting because non-interacting par- 
ticles are always localized in the region of space which 
corresponds to the global minimum of the disorder po- 
tential, i.e. the limit of weak disorder and weak inter- 
actions is already singular. The study of the so called 
'dirty boson problem' was first prompted by the disap- 
pearance of superfluidity of ''He in porous media (vycor) 
[3] . Other condensed matter systems of interest include 
thin superconducting films ;4 , Josephson junction arrays 
[5], superfluid helium films on substrates [S], etc. 

More recently, the realization of ultracold atoms in op- 
tical lattices [3 [5] paved the road to studies of strongly 
correlated systems in a controlled way, and to using them 
as emulators for condensed matter models. In the pres- 
ence of disorder, most cold atomic systems are described 
by the Bose-Hubbard model with disorder in the on-site 
potential. 

Control over disorder has been implemented in optical 
lattices in several ways such as using bicromatic lattices, 
a laser speckle field, and by loading a second (heavy) 
component into the system [9HT3|. Bicromatic lattices 
are produced by superimposing two optical lattices with 
the second one being weak and incommensurate with the 
first one. This method has been used to study Anderson- 
type localization phenomena [10]. Both weakly interact- 
ing and strongly interacting limits have been studied in 
a speckle potential which is implemented by scattering a 
coherent laser light from a rough surface [HI US] • 

Considerable amount of theoretical effort has been ded- 
icated in the past to understanding the phases and phase 
transitions in the system [T^l - Bl] . In addition to the Mott 
insulator (MI) and superfluid (SF) states present in the 
pure case, in the presence of disorder, the ground state 
features a third, non-conducting but compressible phase, 
the Bose glass (BG). Whether there exist a direct transi- 
tion between the MI and SF states has been the subject 
of a long debate. Fisher et al. [2] argued that such a 



transition is unlikely, i.e. for finite disorder MI and SF 
phases are always separated by BG, but alternative pos- 
sibilities were not ruled out rigorously. This resulted in a 
controversy since, on one hand, mean-field type theories 
are inadequate in capturing the physics of rare statis- 
tical fluctuations driving the MI-BG transition, and, on 
the other hand, various numerical techniques are severely 
limited by finite size effects (THl HH H^. The contro- 
versy has been resolved in Refs. [251 [26] which proved 
the theorem of inclusions and concluded that all transi- 
tions between the fully gapped and gapless ground states 
in disordered systems are of the Griffiths type and thus 
the resulting gapless phase is insulating. Moreover, the 
original conjecture OE] that the MI-BG boundary cor- 
responds to the disorder bound, A, equal to the MI gap 
in a pure system turns out to be a rigorous result by 
the same theorem. Here it is important that disorder is 
bounded since otherwise the MI phase is eliminated al- 
together. [We note that proofs based on rare statistical 
fiuctuations are valid only in the thermodynamic limit; 
in finite experimental systems phase boundaries are re- 
placed by crossovers, including complete elimination of 
the BG state for weak enough disorder.] 

In the present work we study a two-dimensional disor- 
dered Bose-Hubbard model and present the first accurate 
results for its ground state phase diagram at unity filling 
factor [27]. The numerical method of solution is based 
on the lattice path integral Monte Carlo using Worm al- 
gorithm [55]. We pay special attention to the critical 
behavior of the system in the limit of weak disorder in 
proximity to the SF-MI transition in the homogeneous 
system. 

The disordered Bose-Hubbard Hamiltonian reads as: 
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where al{a^) is the boson creation(annihilation) opera- 
tor, (ij) denotes nearest neighbor sites, and n — a\ai is 
the density operator. In what follows we use the hopping 




FIG. 1: (Color online) Finite-size scaling of average wind- 
ing numbers squared as a function of disorder bound A for 
interaction strength U /t — 22 and dynamical exponent z = 2. 




FIG. 2: (Color online). Zero temperature phase diagram of 
the two-dimensional Bose-Hubbard model at filling factor n = 
1. The MI-BG transition at A = Eg/2 is obtained using gap 
data from Ref. [31]. (The green triangle is the point on the 
SF-BG boundary obtained in Ref. [28].) 



matrix element i as a unit of energy, U is the on-site 
repulsion between atoms, and /i is the chemical poten- 
tial. Random disorder, e^, is uniformly distributed on the 
[—A, A] interval and no correlations between the sites of 
a square lattice. In our numerical study, we work at fill- 
ing factor n = 1 and, for definiteness, choose the chemical 
potential in the middle of the MI gap when appropriate, 
i.e. gaps for creating particle and hole excitation in the 
MI phase are both equal to Eg/2- 

To construct the phase diagram we employ standard 
procedures. The SF-BG transition lines are determined 
from finite size scaling of the superfluid stiffness calcu- 
lates from the statistics of winding numbers squared, 
(W^), using formula A, = {W^)L^-'^/Td [30]. Accord- 
ing to the theorem of inclusions, the BG-MI line is de- 
termined by the A — Eg/2 criterion (recall, that it is 
impossible to detect this boundary directly in finite-size 
simulations). 

Figure [l] shows {W'^)/2 vs A/t curves at U/t ^ 22 
for different system sizes. In this case we scale space 
and imaginary time dimensions according to the dynam- 
ical critical exponent z — 2 predicted in Ref. [2] (strictly 
speaking, in the thermodynamic limit any value of 2; > 
can be used for determination of the critical point). 
Barely measurable flow of intersection points with the 
system size allows us to estimate transition points with 
relatively high accuracy. From the data in Fig. [l] we de- 
duce Ac/t = 7.76 ± 0.04. The full ground state phase 
diagram in the {U, A) plane is shown in Fig. [2J The 
boundary between the MI and BG phases at A/t = Eg/2 
was constructed using data for the MI gap in a pure sys- 
tem calculated in Ref. [31] . 

Reentrant behavior, similar to the one observed in one- 
and three-dimensional phase diagrams, is also present in 
two dimensions, in agreement with earlier observations 
[161 128|. In compliance with the theorem of inclusions. 



BG always separates SF and MI states which meet only 
at the MI-SF transition point of the clean system at 
Uc/t = 16.7424(5). As we move away in the vertical 
direction, weak disorder works in favor of the SF phase 
shifting the transition points to the right — this appears 
to be a common behavior in all physical dimensions. The 
mechanism, however, is not universal. In ID, it can be 
explained by the destructive interference of the vortex 
instanton contributions to the partition function |20j . 

In the weak interaction limit, C/ — >■ 0, the transition 
line does to zero with an infinite slope, Ac oc vC/ ^32J . In 
this region, interaction is very efhcient in screening deep 
potential wells and stabilizing superfluidity. Numerically, 
this region is extremely hard to study [BSj and we were 
not able to see its asymptotic laws. 

Remarkably, superfluidity persists up to extremely 
large values of disorder and interactions. For interme- 
diate interactions, superfluidity survives when disorder 
potential is about ten times larger than the bandwidth, 
A/t = 70. Likewise, the superfluid 'finger' at A/t = 30 
extends all the way to U/t ~ 48. However, super- 
fluid properties of the system in the large disorder and 
large interaction limits are not robust. In Figs. [3] and 
|4] we plot the superfluid stiffness calculated along two 
representative cuts of the phase diagram, one at fixed 
interaction U/t = 26 and the other at fixed disorder 
A/t = 35. Small values of As for A > 30i ensure small 
superfluid-normal transition temperatures according to 
the Nelson-Kosterlitz relation Tc/t = {tt/2)As{Tc). Note 
that Figs. [3] and |4] can be used to determine upper bounds 
on Tc because As(Tc) < As(0). For example, numerical 
simulations of the SF-normal transition temperature at 
U/t = 26 and A/t = 35 yielded T^/t = 0.07± 0.012 below 
the Tc/t = 0.08 estimate. 

Low transition temperatures have important conse- 
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FIG. 3; Superfluid stiffness as a function of disorder bound 
A at fixed interaction strength, U /t — 26 for system size 
L = 12 X 12. 
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FIG. 4: Superfluid stiffness as a function of interaction 
strength U/t at fixed disorder bound, A/i = 35 for system 
size L = 12 X 12. 



quences for current cold-atom experiments. Though the 
{U/t = 26, A/t — 35) point is chosen to be far from the 
edges of the SF-BG boundaries in Fig. [2] the value of Tc 
is well below typical experimental temperatures which 
are still at (or above) the tunneling amplitude t. Thus 
observing the ground state phase diagram remains a chal- 
lenging task. In current experimental setups, only a small 
fraction of the superfluid region will survive finite tem- 
perature effects. 

Let us now focus on the weak disorder case at the tip of 
the Mott lobe in the pure system. Consider the Euclidian 
superfluid hydrodynamic action [h = 1) 



S ^ dr dTliin)'^ 
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where k, is the compressibility, r is the imaginary time, 
^(r, r) is the superfluid phase field, and {n(f)) is the av- 
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FIG. 5: (Color online). Phase diagram in the vicinity of the 
Mott lobe for weak disorder. MI-BG transition curve is de- 
termined by the MI gap, A = Eg/2, calculated in Ref. [31j . 
Note extremely weak dependence of the SF-BG critical point 
on disorder for small A. 



erage density at the lattice point f (the integral over dr is 
understood as a discrete lattice sum, and the gradient as 
a finite difference) . Since the first term contains full time 
derivative of the phase variable, we observe that (i) the 
integral J $ dr can be only a multiple of 27r, (ii) in a pure 
system with {n(f)) = 1 this term is irrelevant and the ef- 
fective action (after rescaling of the time variable with 
the sound velocity) is that of the classical 3D XY-model, 
(iii) in a disordered system this term is of any importance 
(not a multiple of 27r) only in the presence of topologi- 
cal defects in the phase field. More precisely, its value is 
purely imaginary and for the (2-|-l)-dimensional space- 
time vortex-loop with projected (on the space plane) d- 
dimensional algebraic area A is given by a simple formula 



^loop 



27rj / dr {5n{r)) 



(3) 



where 5n{r) = n{r) — 1 is the local density fluctuation 
about the mean value. 

The nature of the SF-MI transition in a clean system 
described by the 3D XY-model is linked to the prolifer- 
ation of large vortex-loop instantons which disorder the 
phase field. In the absence/irelevance of the phase term 
all vortex-loop instantons act in 'unison', in other words 
their contributions interfere constructively in the parti- 
tion function. This is no longer the case in the disordered 
system since an area integral in Eq. ([3]) is now a ran- 
dom variable. When random vortex phases become large, 
their contributions to the partition function cancel each 
other making them inefficient in destroying phase coher- 
ence across the system. Let us estimate the typical phase 
of a vortex- loop of size ^ where ^ is the correlation length. 
It is proportional to the total particle number fiuctuation 
in the area A ^ ^"^ va response to the random chemical 



potential fluctuation in this region, S^ia = A~^ J CrCpr, challenge, 
which, by the central limit theorem, scales as A/^. Using 
system's compressibility, k — dn/d^, we then find 



Im Sloop ~ 2nAKA/^ -- 2tt^kA 



(4) 



As long as disorder remains perturbative (i.e. S'loop ^ 1), 
we have ^ ~ (1 - UciA)/U)-'' where v = 0.67155 is 
the correlation length exponent and Uc{A) is the critical 
interaction for a given A. Equation Q extrapolated to 
Im S'loop ^ 1 predicts how close U should be to Uc{A) 
to start seeing relevance of disorder, provided the system 
size is exponentially large: L > ^ ex exp(C//A), as it 
follows from the estimate 



Im Sloop (x(A/C/) In e 



(5) 



justified below. 

One might think that Eq. Q implies linear scaling of 
phase with the correlation length. However, in the vicin- 
ity of the particle-hole symmetric SF-MI transition point 
Uc/t, the renormalized compressibility itself vanishes as 
K ^ l/US,, canceling dependence on ^ in Eq. Q. It means 
that each length-scale contributes equally to the vortex- 
loop phase, i.e. the final dependence is only logarithmic, 
bringing us to Eq. ([5|. 

On the other hand, we conclude that d = 2 is the 'crit- 
ical dimension' starting from which the relevance of weak 
disorder can not be seen by assuming that compressibil- 
ity follows the Josephson relation, k cx ^^~'*, meaning 
that the shape of the critical line Uc{A) is associated 
with essentially non-universal microscopic physics. [Fi- 
nite compressibility in Eq. (HI) clearly implies relevance 
of disorder on length scales > (kA)"^/"*.] This is to be 
contrasted to d = 1 where the scale of relevance of weak 
disorder defines the form of the line Uc{A), see pU] . 

In Fig. [5] we show the phase diagram close to the tip 
of the Mott lobe. Critical points were determined using 
finite-size scaling with z = 1. Clearly, critical points for 
the SF-BG transition are at disorder values much larger 
than the Eg/2 boundary for the MI phase. Surprisingly 
enough, the first data points for the SF-BG line are indis- 
tinguishable from the critical value in the clean system 
Uc/t = 16.7424(5) within the error bars. This was unfor- 
tunate to some degree because even though critical points 
were determined with accuracy better than four digits we 
still did not have enough parameter range to make an un- 
ambiguous case for the form of the line Uc{A). 

Summarizing, we have presented the full ground state 
phase diagram of the disordered Bose-Hubbard model at 
unity filling factor. Interestingly, while the superfluid 
phase is remarkably stable against strong interactions 
and disorder it is rather fragile with regards to finite tem- 
perature effects. Our numerical data feature an essen- 
tially vertical SF-BG line for weak disorder; understand- 
ing physics behind this phenomenon in rf > 2 remains a 
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